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ABSTRACT 


t ' ^ 

An extensive 'investigation of the orientation of crystallites in 

"i , * 

molded, artificial graphites has “been' performed. It has "been found 

' ( v - 

that th'e crystallites i-n molded graphites are usually - hut not always - 
distributed symmetrically about some axis. The symmetry axis has been 
found to differ frequently from the molding axis, contrary to what has 
generally been assumed. Furthermore, the direction of the symmetry 
axis and the degree of orientation have been found to vary somewhat 
within a .given graphite block. 

An equation, originally proposed by Pappis, et al. , has been found 
to be suitable for describing the distribution of crystallites about 
the symmetry axis. The equation involves two unspecified parameters. 
Specification of the direction of the symmetry axis involves two 
additional parameters. A technique for evaluating these parameters 
from experimental data obtained by the Bacon method is presented. 

Also presented are recommendations for future research. 
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THE ORIENTATION OF CRYSTALLITES IN MOLDED GRAPHITES 



INTROHJCTION 


Virtually all nonpyrolytic artificial graphites are polycrystalline 

* , r 

and most have their crystallites preferentially oriented in some fashion. 
Because many of the properties of .graphite crystallites are anisotropic 
( i . e ., directionally dependent), preferential .crystallite orientation 
results in anisotropy of the properties of most hulk graphites. The 

t. 

degree of anisotropy of a graphite determines its suitability for' many 

i > , a - * " 

applications. For example, a spacecraft heat shield,, ideally, should 
have a high termal conductivity parallel to the surface in order to 
distribute the heat load over its .entire’ area, but a low thermal 
conductivity perpendicular to the surface to insulate the spacecraft. 
Thus, a graphite intended for use as a heat shield should be highly 
anisotropic. On the other hand, .such a graphite might be unsuitable as 
a moderator in a nuclear reactor because of mechanical problems associ- 
ated with non-uniform expansion due to neutron irradiation. The' degree 
of anisotropy of .a given property is related to the degree' of 'crystal- 
lite orientation of the graphite (see ref. 2) and, consequently, the 
degree of orientation is, itself, an important property. Unfortunately, 
no .simple method for describing, completely and unambiguously, the degree 
of orientation of a graphite is presently available. The purpose of 
this study was to develop such a method. 

Crystallite orientation in artificial graphites results from an 
interaction between the structure of the crystallites and the methods 
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by which graphite bodies are generally manufactured . An extensive 

i / ' i 

, . ] 1 

treatment of graphite structure and. manufacture is- given by Nightingale. 

Q 

Harris, Miller, and Craik discuss the effect of these factors on 
crystallite orientation. A brief discussion of these factors will be 
given here . 

Artificial graphites are generally manufactured from a mixture of 

a granulated, carbonized filler - such as petroleum coke or lampblack - 

and a viscous binder with a high carbon-to -hydrogen ratio - such as 

coal tar pitch. The binder-filler mixture is formed under pressure 

into billets either by molding into a form or by extrusion through a 

die. The raw' billets are then converted to graphite by heating. 

Each grain of filler material contains many regions of partial 

atomic ordering - incipient crystallites - which, after graphitization, 

become .tine crystallites. ^ A graphite crystallite is made up of layers, 

or "basal planes," each of which •constitutes, in essence, a single 

19 

large aromatic molecule. (See fig. 1.) Bonding within each basal 
plane is covalent with a bond energy of 150 kcal per mole; bonding 

between- basal planes i-s due to van der Waals forces and the bond energy • 

11 

is only 1.3 kcal per mole. It is obviously much easier to break bonds 

between basal planes than within them. Consequently, when filler 

material is granulated, the grains tend to form in elongated shapes 

with their major faces parallel to the basal planes of their surface 
9 

crystallites. In general, the smaller a grain is the more likely it 
is that most of its crystallites will have similar orientations, 




although the degree to which this is true depends upon the nature and 

9 ’ 

source of the filler material. 

When the filler -hinder mixture is formed into 'billets under 
pressure, the oblong filler grains experience torques which tend to 
orient them preferentially with respect to the direction of applied 
force. When the forming process is molding, the most stable orienta- 
tion of the grains - and, hence, of the basal planes of their consti- 
tuent crystallites - is perpendicular 1 ' to' the molding direction; in the 

■case of -extrusion, the preferred orientation is parallel to the 

11 

extrusion direction. 

The alignment of the crystallites within each 'grain" and of the 

grains within a graphite body is never perfect, , and thus, although one 

’ *> ' 1 

direction is preferred^ "many crystallites may' be found -in other 

directions. 8. Therefore, a complete and unambiguous description of the 

degree of orientation of a given graphite must include the specif ica- 

b 1 

tion of the relative number .. of crystallites in every direction. 

1 * 

Crystallite orientation is generally determined experimentally by 
means of X-ray diffraction. The principle underlying all X-ray studies 


a The orientation of a crystallite is usually expressed in terms of' 
an imaginary ray normal to the basal planes of the crystallite rather 
than in terms' of the basal planes themselves, as was done above. This 
convention will be used henceforth throughout this thesis. 

^The relative number of crystallites in any direction is the 
actual number in that direction divided by the number in some reference 
direction. 




of orientation is that the intensity of diffraction from the "basal 

Q • ' 

planes oriented in a given direction is proportional to the number 
of' such planes and,, therefore, to the relative number of crystallites 
in that direction. A' number of X-ray techniques have "been pro- 

*1 _ II *7 1 n 1 l y , f 

posed. ’ ’ The principal technique utilized In this study 

» 

p 

is that of Bacon, but some data were obtained by .the method of AH, 

\ ' 1 

Fitzer, and Ragoss.^ 

For convenience, this study was limited to molded' graphites, but 

► , / , 

of describing orientation which is developed should be 

- # 

to extruded graphites as Well. 


the method 
applicable 


a The basal planes of graphite are frequently designated in X-ray 
diffraction work as (002) planes, the numbers in parentheses being 
Miller indices. Both designations are used interchangeably in this 
thesis. The diffraction of X-rays from the basal planes is commonly 
referred to (002) reflection, although this use of the term "reflection” 
is not rigorously correct-. 



SURVEY OF PREVIOUS WORK 


The earliest systematic study of the orientation of crystallites 
in graphite was reported by G. E. Bacon. 2 The experimental arrangement 
for the method of Bacon is illustrated in figure 2. The incident X-ray 
beam, x, is horizontal; the X-rays are unfiltered CuK radiation. The 
graphite specimens are thin flat plates approximately 2 X 1 X 0.1 cm, 
cut with their 1 cm edges parallel to the molding axis, P. Each 
specimen is aligned so that its 2 cm edges are horizontal and perpen- 
dicular to the incident X-ray beam, and its 1 cm edges are tilted from 
the vertical,- z, by 13° * a 

The graphite specimens are thin enough to allow transmission of a 
significant fraction of both the diffracted and undiffracted components 
of the X-ray beam. The diffracted X-rays form a cone with a half -angle 
of 26° which is twice the Bragg angle for the (002) reflection using CuK 
radiation. When the diffracted X-rays strike the photographic film - 
which is perpendicular to the undiffracted X-ray beam and, thus, to the 
axis of the cone - they form a circular image of varying density on the 
film. The density of the image at any angle, i, is proportional to the 
number of crystallites at some related orientation within the specimen. 
The film density is determined at each angle of interest with a micro- 
densitometer . 

a This tilt is necessary so that crystallites with orientations 'Of 
0° to 13°, with respect to the forming, axis, can be detected. . A 
theoretical discussion of this point is presented in the section, 
"Theoretical Analysis." 
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The main theoretical' problem of tlie Bacon method is the determina- 
tion of the' relationship between angles on the film image and the ■ 
orientation of crystallites within the graphite specimen. Bacon 
assumed that the crystallites in artificial graphites are symmetrically 
distributed about the forming axis. Such a distribution would mean 
that the only parameter essential to. a description of the orientation 
of the crystallites is -the angle, 0, which their normals make with 
the forming, axis. Based on this assumption, Bacon derived the follow- 
ing equation relating 0 and § : 

cos 0 = cos 2 77° - sin 2 77 9 sin £ (l) 

For convenience in expressing the angular distribution of crystal- 
lites,. Bacon defined an orientation function l(0) equal, to the 
relative number of crystallites per unit solid angle about the inclina- 
tion angle 0. l(0) is usually normalized so that 1.(0’) = 1 when 

0=0. In practice, l(0) is taken as the relative diffraction intensity 
which, in the case of the Bacon method, is assumed to be equal to 
D(0)-/d(O) , where D'(0) and B(.0) are the' film densities at 0 and 0°, 
respectively. 

The Bacon method of determining crystallite orientation has the 
advantage that with one specimen and one exposure it furnishes- a con- 
tinuous, quantitative mapping of l(0) . The principal disadvantage of 
the method is that it does not yield diffraction intensities directly 
and immediately, but requires the intermediate steps of film processing 
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and microdensitometry. These steps are not only inconvenient , but also 
constitute possible sources of error. 

Ali, Fitzer, and Ragoss^ contend that a plot of l(</) versus (j> in 
polar coordinates will always be elliptical within experimental error. 
Based on this contention, they have proposed that intensity measurements 
need be made at only two angles, 0° and 90°, since' the value of l(0) 
for any intermediate angle can be calculated from the polar equation 
jf an ellipse 


■I($ = 


1(0)1(90) 

[l(0) 2 sin + l(90) 2 cos 2 ??] 1 / 2 


( 2 ) 


The experimental procedure proposed by Ali, et al., is as follows: 
Two flat plates, or discs, are cut so that one is parallel and one is 
perpendicular to the forming axis of 1 the graphite to be "studied. Each 
specimen is mounted in the flat specimen holder of an X-ray^ diffracto- 
meter and, using a counter-goniometer, the (002)’ peak of each is 
determined and recorded on a strip chart ‘recorder-. .The value, of l(90) 

, V ♦ ‘ . 

is taken as the peak height of t the ' specimen whose face is parallel to 

the forming axis (0 = 90°) divided by' the peak height of the perpendicu- 

1 ‘ * % 

lar specimen = 0°).. The value of l(0) . is taken as one. If 
intermediate angles are to be investigated, a separate specimen must! 
cut for each. 

The method of Ali,. et al., is quite simple "if only two specimens 
need be cut and studied. However, this will be the case only if the 
distribution of crystallites in the graphite to be studied is, in fact, 



II 


"both symmetrical and elliptical about the forming axis.* 'Unfortunately, 

« » > * ' » 

Ali, et al. , present only limited experimental substantiation of their 

* 

assumption and they give no quantitative results. , 

Karris,, Miller, and Craik^ present polar ( plots of 1(0) versus (j 
for four graphites, and in no case is the distribution elliptical. How- 
ever, all four graphites were extruded rather than molded and three of 
them were specially prepared in the laboratory rather than commercially 
manufactured. These results, therefore, do not rule out the possibility 
that molded commercial graphites possess elliptical distributions. It 
should be noted that if All's assumption of a symmetrical, elliptical 
distribution is correct, then the normalized value of 1(90.) constitutes 
a single parameter which, in conjunction with equation (2), provides a 
complete specification of' the spatial distribution of graphite crystal- 
lites. 

An alternative one-parameter equation that has been proposed by 

278 i8 a 

several investigators is ' ' ’ ’ 

1(0) = cosM0 (3) 

This - equation works fairly well with pyrolytic graphites, but it 
usually does not work well with molded graphites because they generally 
have some crystallites with their normals at 90° to the molding axis. 
Equation (3) assumes that l(90) = 0. 

a Actually, Bacon proposed the form l(0) = sin^0 which was intended 
to be applicable to extruded graphites. The direction of maximum 
orientation of molded and extruded graphites differs by 90° • 
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Pappis, et al., 1 ^ have proposed a variation of equation (3) thal 
overcomes this problem. 

1(0) = A eos^ + B (*0 

Unfortunately,- these investigators also fail to present any quantitative' 
substantiation of their equation with experimental data. The equation 
is still of considerable interest, however, and should be compared with 
the ellipse proposed by Ali, et al., to see if either is clearly superior 
to the other. Equation (4) can be regarded as a two-parameter equation 
since, by proper normalization of l('0)>, we may let A + B = 1. 

Several investigators 1 ' 2 -> 10 12 have proposed orientation param- 
eters which attempt to specify the degree of orientation of a graphite 

with a single number. Such numbers, although useful for some purposes,, 

\ 

do not describe the spatial distribution ’Of : crystallites. Therefore, 

they will not be discussed further here.’* 

The assumption of a symmetrical distribution of crystallite's about 

the forming axis has been mentioned repeatedly in this section and in 

* • 1 v . 

the Introduction, and it is either seated or implied ( in .most papers on 

' ' l . ■ * 

graphite orientation. ^.>2, 9,10, 12,13 >15.> 16 Cavin, 10 however, has 

" > x ' ) < " i * 1 

‘recently reported experimental results which' contradict this assumption. 
In his investigation, Cavin observed -a shift' of the symmetry axis from 
the forming axis by as much as 12°.' The experimental technique used by 
Cavin was a form of a pole a figure technique employing a Schultz 

a A "pole" is an imaginary normal to a crystallographic plane; a 
(002) pole is thus identical to the basal plane "normal" utilized 
throughout this thesis. 
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preferred-orientation apparatus . ^-7 The "basis' of this, and all pole 
figure techniques, is the determination of the diffraction intensity not 
only as a function of the angle of inclination,, with regard to 
some reference axis hut also' as a function of the aximuthal angle, 
t), lying in a plane perpendicular to the reference axis . ^ The methods 
hy which pole figures are obtained and interpreted are considerably 
more complex than the methods proposed by Bacon and Ali, et al. 
Nevertheless, if the .methods of Bacon and All are inadequate to 
determine completely crystallite distributions, the use of a pole 
figure technique may be necessary. One purpose of this study was to 
determine -whether the Bacon or Ali techniques can be modified so as 
to detect and account for possible -shifts in the symmetry axis and, 
thus, to yield a complete description of the distribution of. crystal- 
lites in molded graphites. 


a ln Gavin's study, the reference axis was the forming axis. 



EXPERIMENTAL TECHNIQUES 


p 1 

The methods of Bacon and All, Fitzer, and Eagoss were employed 
exclusively in this study. The vast majority of the data were taken 
using the Bacon method, hut a few checks were made with the Ali method 
becaus'e of its simplicity. Each -of these methods is outlined in the 
previous section. This section will present details of the methods as 
employed In this study. 

The X-ray instrument utilized for both the Bacon and Ali methods 
was a General Electric Co. XRD-5 with various accessories. The X-rays 
were nickel-filtered CuKa radiation. The peak tube voltage was 50 kV. 

The experimental apparatus used for' the Bacon method was a trans- 

’ « * * 

. . * ' 

mission Laue camera with 0.020 -inch-diameter pinhole collimator and 

a 4- by 5-inch film cassette-. t Ten grades of graphite were investigated 

' • . . ■ - . a 

by the Bacon method. Three specimens of each grade were studied. 

r 4. , f 

d , * 

The grades studied are listed in table 1 along with certain of their 
properties. 

The specimens, which were 2-5 X 1 X 0.1 cm, were mounted in a 
specially constructed holder which permitted them to be oscillated 
horizontally in a plane normal to the incident X-ray beam. The 
oscillation increased the number of grains irradiated. The period of 
oscillation was 1 minute and the amplitude was 0.75 inch. The speed 

Only two specimens each of grades ATJ and 2 d8d were available. 
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TABLE 1 


SOME PROPERTIES OF GRAPHITE GRADES STUDIED 


Grade 


■ Density , gm cm"^ 

' i 

Maximum grain 
size, mm 

ATI 

n 

PC b 

1.74 

O'. 15 

ATJ (G.P.) a 

PC 

.1.72 

0.15- 

2BE 

PC 

l.,4o 

. 0 . 15 ' 

9RL 

PC 

1.68 

,o.oer 

3^993 

PC 

1.63 ‘ 

0.08 

•1(007 

PC 

1-70 

0.20- 

L31 

LB C 

1.66 

0.15 

2D8D 

LB 

1.1(0 

0.18 

GDG 

PC and LB 

.1.1(9 

^ 0« ^KL 

CMB • 

PC and LB ‘ 

1.79 

0.08 


Gas purified. 
^Petroleum coke. 
c Lamp'black. 






of oscillation was kept constant so that all grains irradiated were 
exposed for' the same length of time. - The exposure time was 1 hour. 
Specimens of grade CDG, which is somewhat coarse grained, were exposed 
in two steps of l/2 hour each. _ After the first l*/2 hour, the specimens 

f ' i 

were raised in the holder so that more grains •would he exposed. 

The film used was Kodak Type M industrial X-ray film. The films 

were individually processed with Kodak Liquid X-ray Developer and 

* 

Replenishes and Kodak Liquid X-ray Fixer and Replenisher. The 
manufacturer 1 s processing instructions were followed throughout. 

The density of the processed film was determined with a Joyce Loehl 
dual -beam microdensitometer-. The films were mounted on a special 
rotary stage with angular markings to 0.1°. The stage was then rotated 
to each desired value of £ and radial scans of the circular diffrac- 
tion image were made. The resuLt of each such scan was a trace, -on 
ruled paper, of the density versus tan 20. Typical traces are shown 
in figure 3- 

Bacon stated in his paper that the integrated density - the area 
under the diffraction peaks - is proportional to !(•$) . All, et al. , 
and Guentert contend that the peak height and peak area are proportional 
to each other and to l($) and that the height is preferable to the 
area since it is easier to measure. Actually, the height and area are 
not always proportional to each other as is shown in table 2 in which 
the ratio of the normalized peak heights and peak areas is listed for 

a 9 is defined as the angle between the incident X-ray beam and the 
basal planes of the diffracting crystallites. The angle between the 
diffracted and undiffracted components of the X-rays is 29. 














the two specimens of grade ATJ. It can he seen that the ratio is not 
unity and,, in fact, is not even constant. The variations are too 
large and systemmatic to he due entirely to experimental error. 

The problem of ascertaining the best measure of l(0) is compli- 
cated by the following facts: (l) The size and shape of a diffraction 

peak are influenced by other factors in addition to the number of 
crystallites causing diffraction (see, for example, ref. 5)* (2) A 

molded artificial graphite is not a homogeneous material. It has at 
least two phases - binder and filler - and more than two phases if two 
or more fillers are used. The overall diffraction peaks obtained are 
the sums of the peaks of the individual phases. This problem is treated 

1 p 

in some detail by Noda and Inagaki. (3) The abscissa of the 
diffraction peaks obtained from Bacon films is tan 20 rather than 
29. 

The resolution of this problem was beyond the scope of this study. 
Furthermore, it was felt that a comparison of the results obtained 
using peak height and peak area data would be of interest. Therefore, 
both were measured in this investigation.. 

The 28 films obtained from the 10 grades of .graphite studied were 
scanned in 10° increments of £ from 0° through 350 °; In addition, 
one grade - 2BE - was reread every 2.5° of £. The heights of the 
diffraction curves were obtained by subtracting the average* background 

„ k 

height from the average height of the peak crests, both of which were 

i . 

read directly from the ruled paper. 



Three grades of graphite - ATJ, 2D8 d, and GMB - were investigated 
"by the Ali technique. The ‘-approximate dimensions of the specimens used 
were 2-3A X 1 X 5/16 inch. For each specimen, the value of 28 was 
varied continuously over about a 5 range encompassing the Bragg value. 
■The variation was extensive enough to encompass the entire peak as well 
as some background on either side. The peak area was determined with a 
planimeter and the height was taken as the peak deflection minus the 
average of the background deflections. 

The peak heights measured in this, study have .an uncertainty of 
±1 to 2 percent* 3 except when a dust speck or film blemish caused a 
spurious deflection at or near the peak crest. In such cases the 
uncertainty is est ima ted to be about ±5 percent. . The uncertainty in 
peak areas is considered to be .greater than in peak heights because of 
the increased effect of uncertainties in the baseline. Errors in 
determining the baseline have a linear effect on the uncertainty of the 
peak heights but a much greater effect on the uncertainty of the peak 
areas because of the divergence of the diffraction peak at its base 
(see fig. 3) . Dust specks and film blemishes are also troublesome when 
areas are being determined. All such spurious deflections were faired 


a The goniometer dial of a diffractometer generally reads values of 
29 rather than 9. The Bragg value of 20 for graphite with CuKd 
radiation is about 26 - 5 °; the value varies somewhat with crystallite 
size and degree of graphitization . 

^The background and peak crests could each be measured to within 
about half of a line on the ruled paper. Most of the peak heights 
obtained lie in the range 50 to 100' lines. 



through before the areas were determined, hut some additional uncertainty 
was introduced. Overall, the planimeter integrations are' estimated to 
he uncertain to about 2 to J percent. 

Many factors contribute to the uncertainty in the angular orienta- 
tion of the specimens, but the major source of uncertainty unquestionably 
was the cutting process. The Bacon specimens were cut. in two steps: 
.First, 1-inch cubes were sawed from the billets as received at this 
laboratory. The cubes were cut with an angular accuracy of about 2° 
or better and the pressing direction was clearly marked. The specimens 
were then cut from the cubes by the Speer Carbon Company under . contract 
to NASA. No estimate of angular uncertainty was furnished by Speer. A 
value of 2 percent will be assumed since the final cutting of the speci- 
mens should not have been any more inaccurate than the sawing of the 
cubes. All other sources of error contributed less than 1°. Thus, the 
total angular error in the Bacon specimens should be less than 5 °- The 
Ali specimens were cut in one step at this laboratory; their total 
angular uncertainty is no more than 2°. 



THEORETICAL ANALYSIS - 


In this section we will consider the- implications' of assuming that 

» * » 

the distribution of crystallites in molded graphites is symmetrical 
about some axis which is not necessarily coincident with the molding 
axis. We will then attempt to devise a technique by which (l) the 
spatial distribution of crystallites in a given graphite can be 
completely determined from a properly obtained Bacon film image, and 
(2) the distribution .so determined can be expressed, completely and 
unambiguously, with a minimum number of parameters. 

Assume that the graphite specimen shown in figure 2 possesses a 

A 

symmetry axis having some unspecified direction., Let S be a unit 
vector coincident with the symmetry axis. In the coordinate system 
shown in figure 2 

A A A A 

S = cos qgi + cos Pgj + cos 7 gk 

= 2 s i + mgj + n s k (5) 

where ccg, 0g, and 7g are the angles which the symmetry axis makes 

, and n^ are the 

corresponding direction cosines. 

A 

Let W be a unit vector normal to the basal planes of some crystal- 
lite of interest’. In terms of the direction angles '( 3 ^, and 7 ^ 

and the direction cosines l , m, and n 

w w n 


with the x, y, and z coordinates axes and l, 


S’ 


“s 
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/N A ’ J V A] 

N = cos ct^-i + cos (B^j + cos y-^k 


= V * V + “H* 

(6) 

/s 

The scalar product of the unit vectors N and S is, 
tion, equal to the cosine of the angle between them. But 

A A 

between if and S is the orientation angle (fi . Thus, 

by defini- 
the angle 

A ,A 

N • S = ZjjZg + n^mg + n N ng 
= cos 4 

(?) 

From the properties of direction cosines, we know that 


2 2 2 
*** * n N “ 

(8) 

p p p 

V * *s + n s ’ 1 

(9) 


If we solve equation (8) for mjj and equation (9) for ng, we get 




We may adopt the convention that, the vector S always has an upward 
component so that equation (ll) "becomes positive only. 

The condition for (002.) diffraction to occur is that the angle 
between the hasal planes and the incident X-ray. beam' must be the Bragg 
angle for the radiation used. The Bragg angle for (002) diffraction 
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using GuKa, radiation is 13.25°- If the basal planes make an angle of 
13 . 25° with the incident X-ray team, their normal must make an angle 
of either 76 . 75 ° or 103 - 25 °, depending on whether the x-component of 
the normal is positive or negative. If we require that it "be positive, 
then 

ojj- = 76.75° ( 12 ) 

It can "be shown from consideration of spherical coordinates that 

n^ = sin 76 . 75 ° sin | (15) 

If equations (10) through ( 13 ) are combined with equation ( 7 ), the 
following general equation relating 0 and. £ is obtained: 

cos 0 = cos 76 - 75 ° ± rn g sin ' 76 . 75 ° \fl - sin 2 £ 

+ y i - m Q 2 sin 76 . 75 ° sin i W 

Note that equation ; contains vwp uiuwiuwii jjoiajucuciD, Iq and 

* 

mg, the direction cosines of the symmetry, axis with respect to the x- 

' , , , i 

and y-coordinate axes. -Bacon, in, assuming that the symmetry- axis is 


a The second term on the right-hand side of equation ( 1^-) is nega- 
tive when 90° < i < -270°, otherwise it is positive. -This results from 
the fact that for the diffracted beam to strike the left half of the 
film, the y-component of the crystallite "normal vector must be negative. 
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coincident with the forming axis, - in- effect ass-umed that ctg = 105-25 
and p g = 90 °. This is equivalent to the assumption that 
l Q = -cos 76.75 a and mg = O'. If these values of and mg are 

inserted into equation (15 ) > the result is 

cos (f> = -cos 2 77° + sin 2 77° sin £ (15) 


Equation (15) is essentially the equation derived by Bacon; it differs 
from Bacon’s- equation only in the sign of the terms on the -right. Since 
we are only ‘concerned with values of (f) between 0° 'and 90°, we may 
require cos <j to be positive and rewrite equation (15) as 


cos 


4 - : 


cos 2 77° - sin 2 77° sin £ 


(16) 


Equation (l6) results in all values of <J) from 0 ° to 90 being 

represented on the film image as is indicated schematically in figure 

* * (* . ^ *■ 

4(a). Many values of 4- are represented at several points on the 
film. Figure 4(b) is a schematic representation of ,the limiting value- 
of (ft and £ if the specimens are aligned vertically rather than 
tilted by 15°, as specified by Bacon. These values result from- -apply- 
ing the condition that ctg = 90°, and therefore, ig = 0 to equation 
(l4). It is obvious from figure 4(b) that crystallites, with orientation 
■angles between 0° and 15° would not be detected if the specimens were 
aligned vertically. 


^acon took the Bragg angle to be 15°, and thus lg would be 
- cos 77° • 
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If Bacon 1 s assumption of the coincidence of the symmetry and form- 
ing axes .is incorrect, two difficulties arise: (,l) the values of ^g 

and nig are not known a priori , (2) not all inclinations, <j), ar.e ^ 
necessarily represented on a given film image. These difficulties cai 
he overcome, at least in principle, if a valid closed-form expression 

7 , t 

relating I and 4 is known and if sufficient experimental data are 
available. The procedure involved will be illustrated by considering 
equation (4-'), which is l($) = A cos ^4 + B. -Note that (fs appears on. 
through its cosine. Equation ( l4) gives cos 4 as a function of the 

'* * , I 

parameters l Q and m and the variable £ . If equation's (4) and (l4) 

D O' 

are combined, we get 

T = A^g cos 76.75° ± mg Sin 76-75° \/l - sin^fi 

+ ^1 - l s 2 - mg^ sin 76.75° sin + B (17) 

Equation (17) contains five unknown parameters - A, B, M- 

and mg - which can be determined, in principle, if five or more pairs 

* 

of experimental values of -I and £ are available. Appendix A 
describes a least squares method for solving nonlinear equations such 
as equation (17). Appendix B gives a listing of a computer program and 
related subprograms by means of which the solution can be carried out. 
This method and program were used in this investigation. 

a The success of solutions of equation (17) and the subsequently 
derived equation (19) is dependent on the number and accuracy of the 
data points available. 



Equation (2) can also be combined with equation (l^) to yield an 
equation si m ila r to equation (17)* This can be seen more easily by 
rewriting equation (2) as 


I = 



CL 

(18) 


This expression may be inserted into equation (l4~) to yield 


I = 


ab 


[a 2 + (b 2 - a 2 )(2g 2 cos 76.75° 


± mg sin 


76 . 75 ° \j: 


1 - sin 2 ? 


nT 7 


Iq - mg sin 76.75° sin i~}' 


l/2 


( 19 ) 


Equation (19) can be solved by the same least squares method applied to 
equation (.17) to yield the best values of a, b, Zg, and mg. The 
requisite computer programs and subprograms are listed in appendix B. 

Once the parameters in either equation (17) or equation (19) have 
been determined, the value of I at any angle (fi can be calculated 
(using, eq. (1*0 to relate (f> and £)• Thus, it is not necessary that 


The parameters l(o) and l(90) have been replaced by a and b, 
respectively, for the sake of generality. We have no justification for 
attaching more weight to experimental values measured at 0° and 90° than 
to values measured at other ans-les. 



29 


a particular angle "be studied' experimentally, or even that it "be 
represented on the film image. 

The number of parameters actually specified need he only four in 
the case of equation (17) or three- for equation (19)-, since, hy proper 
normalization, A + B = 1 and a = 1. The necessary normalization is 
performed in the- programs listed in appendix B- Equations (17) and (19) 
can he compared hy applying each to several sets of experimental data 
and comparing their variances. 

The angle hy which the symmetry axis is displaced from the molding 
axis can he. obtained hy first representing the molding axis hy the unit 

/s 

vector F where 

/\ o ~ _ _o r 

P = cos 105 i + cos 13 k • 

?= -cos 77° i + cos 13° k (20) 

A A 

and then taking the scalar product of S and P 

^ ^ o / “2 2 -1 -2° 

S • P = -l s cos 77 + V 1 - V - mg cos 13 

' 4 

= COS -5 (£l) 

1 

The arcosine of cos 6 is the desired displacement angle, 5. 



RESULTS AND DISCUSSION 


The data from the 28 Bacon-method diffraction films were reduced 
using equations -(17) and (19). The least-squares values of the various 
parameters contained in these equations were determined for each set of 
data using the technique outlined in appendix A and the computer programs 
and subprograms listed in appendix B. The values of the parameters 
which were computed are presented in table 3 . Also presented in this 
table are the values of 5, the angle between the symmetry and reference 

-u £> 

axes, and cr , the variance of the data.- 

Three important observations can readily be made from the results 
presented in table, if: (l) the values of the parameters obtained using, 

peak heights to represent 1(0) generally do not agree with the values 
based on peak areas, even, within experimental uncertainty. (2). The 
values of the parameters computed using equation (17) frequently do not 
agree with the 'values obtained using equation (19)* (3) The values of 

S, in many cases, are too large to be accounted for by experimental 
error. Each, of these points and others related to them will be discussed 
in this section. 


^he values of the parameters a and A are not presented in 
table I since the data and results were normalized so that a = 1 and 
A = 1 - B. Although equations (17)’ : and (19) were solved in terms of 
the direction cosines' Zg and .mg, the direction angles ctg and 3g 
are easier to visualize-; therefore, i these - angles rather than the 
direction cosines. are presented in table I. 

< ,v ' T ‘ • 

^The’ reference axis was the molding axis within the limits- of 

» lr * - 

experimental error. 



TABLE 3 

LEAST SQUARES VALUES OF ORIENTATION PARAMETERS 
(a) Equation (17) - Peak Area Data 


Grade 


Specimen 

number 


2BE 

2HE 

2BE 

9RL 

9RL 

9RL 

3^99S 

3^99S 

3499S 

4oo7 

4oo7 

4oo7 

L>1 

L31 

L31 

2D8d 

2D8d 

CDG 

CDG 

CDG 

CMB 

CMB 

CMB 


a Gas purified. 
b No solution. 


B 

M 

a s _, aeg 

P s ; <3.eg 

5, deg 

a 2 x 10 4 

0.52 

2.8 

114 

91 

11 

2.1 

0-52 

2.6 

122 

92 

19 

l.l 

0.42 

3.1 

122 

89 

19 

2.2 

0.44 

2.9 

118 

90 

15 

2.7 

0.42 

3.1 

114 

91 

11 

3.6 

0.35 

3-0 

ill 

88 

8 

19.5 

0.31 

3-3 

103 

84 

6 

1.7 

0.33 

3.2 

101 

91 

3 

8.1 

O .52 

2 -7 

106 

91 

3 

2.4 

0.50 

2.6 

120 

89 

17 

2.2 

0.52 

2.5 

122 

87 

19 

2.0 

0.59 

2.3 

117 

87 

15 

2.6 

O .52 

2.5 

108 

88 

5 

2.4 

0.52 

2.6 

111 

93 

8 

2.9 

0.50 

2.9 

110 

100 

12 

. 2.2 

0.51 

2.5 

116 

78 

19 

3-0 

0.43 

2.8 

112 

90 

9 

2.7 

N.S. b 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

0.8l 

2.6 

149 

95 

46 

1.9 

N.S. 

N.S. 

N.S. 

N.S. 

, N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

0.6l 

3.0 

107 

78' 

13 

3.1 

0 . 60 ' 

2.5 

107 

102 

12 

2.7 

0.59 

2.6 

128 

87 

25 

2.6 

0.66 

2.2 

124 

94 

22 

3.6 

0.67 

2.1 

120 

88 

17 

6.3 

0.65 

2.3 

117 

90 

l4 

9.4 





















Grade 

Specimen 

number 

B 

M 

a S: deg 

3s, 3eg 

8, deg 

a 2 x 10^ 

ATJ 

1 

0.48 

2.6 

103 

90 

1 

1.1 

ATJ 

2 

0.51 

2.7 

109 

91 

6 

0.7 

ATJ (G.P.) a 

1 

o.4i 

3-0 

111 

89 

8 

1.8 

ATJ (G.P.) 

2 

0.43 

3-0 

108 

90 

5 

2.6 

ATJ (G.P.) 

3 

0.38 

3-0 

107 

90 

4 

1.6 

2BE 

1 

0.34 

3.0 

104 

88 

■ 2 

21.5 

2BE 

2 

0.29 

3.4 

98 

84 

8 

2.0 

2BE 

3 

0.31 

3-2 

95 

91 

8 - 

7.8 

9EL 

l 

0.48 

2.6 

97 

90 

6 

1.5 

9RL 

2 

0.49 

2.5 

109 

88 

7 

1.1 

9RL 

3 

0.51 

2.6 

108 

87 

7 

1.1 

3499s 

l 

0.58 

2.5 

101 

87 

4 

2.7 

3^99S 

2 

0.50 

2.7 

97 

87 

7 

1-7 

3^99S 

3 

0.48 

2.7 

101 

91 

2 

0.6 

4007 

1 

0.46 

3.0 

102 

98 

8 

1.8 

4oo7 

2 

0.48 

2.6 

106 

T9 

12 

1.3 

4oo7 

3 

0.4-1 

3.1 

101 

89 

3 

1.5 

L31 

l 

0.86 

2.2 

109 

86 

8 

1-3 

131 

2 

0.87 

1.9 

■ 102 

91 

2 

l.l 

L31 

3 

0.86 

1.7 

115 

87 

12 

1.4 

2d8d 

1 

0.83 

1.9 

116 

88 

13 

1.8 

2D8D 

2 

0.86 

1.7 

113 

90 

, 10 

1.7 

CDG 

1 

0.57 

2.5 

99 

80 

11 

1.2 

CDG 

2 

0.57 

2.9 

96 . 

100 

12 

1.1 

CDG 

3 

0.64 

2.3 

110 

86 

9 

3.3 

CMB 

1 

0.67 

2.4 

108 

94 

7 

2. '4 

CMB 

2 

0.63 

2.4 

107 

88 

4 

1.4 

CMB 

3 

0.62 

2.7 

105 

90 

2 

■ 1.8 


a Gas purified 
^No solution. 
























Grade 

Specimen 

number 

b 

a s j. deg 

P s > <aeg 

deg 

cf 2 'X kA 

ATJ 


0.51 

109 

91 

7 

4.2 

ATJ 

2 

N.S. b 

N.S. 

N.S. 

■ N.S. . i 

N.S. 

ATJ (G.P.) a 

1 

jr.s. 

N.S. 

N.S. 

N.S. 

N.S. 

ATJ (G.P.) 

2 

n.s. 

N.S. 

N.S. 

N.S. ' ! 

N.S. 

ATJ (G.P.) 

3 

o.4o 

108 

' 91 

7 

17.1 

25BE 

l 

0.37 

104 

88 

2 

35.9 

2BE 

2 

0.31 

96 

84 

9 

21.6 

2BE 

3 '1 

0.32 

95 

91 

8 

25-1 

9RL 

l 

0.51 

101 

91 

2 . 

6.0 ' 

9RL 

2 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

9RL 

3 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

3^99S 

1 

0.59 

112 

87 

10 

5.4 

5499S 

2 

0.52 

102 

88 

2 

6.8 

3499S 

3 

0.52 

105 

93 

4 

8.3 

4007 

1 

0.50 

105 

100 

11 

6.7 

4007 

2 

0.52 

111 

77 

16 

7.0 

4007 

3 

0.44 

105 

90 

2 

12.1 

L3l 

l 

U.s. 

N.S. 

N.S. 

N.S. 

N.S. 

L3l 

2 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

L31 

3 

N.S. 

N.S. 

, N.S'. 

N.S. 

N.S. 

2d8d 

1 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

2D8d 

2 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

CDG 

1 

0.60 

105 

77 

13 

3.7 

CDG 

2 

0.60 

103 

102 

12 

4.6 

CDG 

3 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 

CMB 

1 

N.S. 

! N.S. 

N.S. 

N.S. 

N.S. 

CMB 

2 

O .67 

; 115 

87 

13 

6.7 

CMB 

3 

N.S. 

N.S. 

N.S. 

N.S. 

N.S. 


a Gas purified 
b No solution. 









































The fact that the peak height and peak area data .generally do not 

yield the same values for the various orientation parameters conflicts 

1 7 

with the assumption of All, et al. , and Guentert that they -are 
equivalent measures of l(0) . Although this investigation did not 
attempt to resolve the, question of which is, theoretically., the "better 
measure of l(0) , the peak height data have "been shown to be more 
consistent. For example, every one of the 56 cases listed in table 3 
which involve peak height data converged to a valid solution, whereas 
17 of the 56 cases involving peak area data failed to converge. Also, 
the average variance of the cases involving peak heights is less "than 
the average for the cases involvi-ng areas by a factor of 0.8l.‘ a 

Similar comparisons can be used to show that equation (17) yields 
more consistent results than does equation (19) • All but four cases 
involving equation (17) converged, but 13 eases involving equation (19‘) 
did not converge. The average variance using equation (17) is only 
0.31 of the average- variance using equation ( 19) ■ Furthermore, the 
individual variance for equation (17) is smaller in every case except 
one, for which it is equal to that obtained with equation (19)* 

These results indicate that the form of equation proposed by Pappis, 
et al.^, is better able to describe the spatial distribution of 
graphite crystallites than is the elliptical form proposed by Ali, 

1 ■ l 

et al.\ 

a The average variances mentioned here and in the next paragraph do 
not include the variances of specimens No. 1 and 2 of grade 2BE-, which ^ 
are believed to be atypical. This point will be discussed later in this 

section. 
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Perhaps the most significant observation to be' made from table 3 
is that in well over half of the- cases the value of 5 is greater than 
5°; which is the estimated maximum angular error. In several cases, 6 
is more than 10°. Thus, significant angular differences between the 
pressing axis and the experimentally determined symmetry axis exist in 
many of the specimens investigated. This corroborates the finding of 
Cavin^ and contradicts the assumption of Bacon 2 and others. ^’9, 10, 12, 13,15,16 
The importance of experimentally determining the symmetry axis and 
using it, rather than the pressing axis, as the reference axis for 
specifying the orientation angle, 0, is illustrated by figure '5- In 
this figure the peak height data of film Bo. 201 are plotted in' the- 
following two ways. In part (a), 0 -is referred to the pressing axis;. ■ 

in part (b), ( is referred to- the computed symmetry axis. It is 

obvious that the data points are much less .scattered when • $ is 
referred to the symmetry axis. 

It is noteworthy that 6 is generally not .constant -among the three 
specimens of each grade of graphite. These variations in 5 are some- 
what confused by the fact that for each specimen as many as four values ' 

i , • 

of 5 were obtained which generally differ among themselves..* ’ Tlie most 

consistent values of the various parameters appear to be those- obtained 

* , ' 

from equation (17) using peak heights. If this set of values of . 8 is 
considered, angular differences of 5° or more exist among the specimens 
of 6 of the 10 grades studied. It appears that the symmetry axis of a 
graphite body does not necessarily have the same direction at all points 
within the body- Furthermore, the degree of orientation was also 
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found- to vary throughout a graphite "body since the parameters b , B, 

<i • ! 

and M generally vary somewhat among the specimens within a given 
grade. It is not known why the direction of. the .symmetry axis and the 
degree of orientation vary within a graphite "body , hut the answer 
probably involves an uneven distribution of forces within the body 
during the forming process. Unfortunately , the positions of the 
specimens within the body were not noted and, thus, no patterns regarding 
the variations can be ascertained. 

Two of the specimens of grade 2BE behaved quite differently from 
the other specimens tested in that the variances associated with them 
were abnormally high. Furthermore, the deviations between the experi- 
mental and calculated values of l(0) (based on the least-squares- 
parameters) were not random. This is shown in table U. The possibility 
was. considered that perhaps the values of the -orientation parameters 
obtained were not the true least-squares values, but rather were 
spurious values to which the solutions of equations (17) an^ (19) bad 
converged. Convergence to unreasonable solutions sometimes occurs 
with the technique used if the initial values of the parameters chosen 
are too far removed from the correct values. To check for this, the 
initial values of the parameters ' were systematically "varied over a 
wide range, but no solution other than that presented in table k was 
obtained. The three films for this grade were then reread at intervals 
of £.5° in i to minimize the effects of errors in the individual data 
points. The data were reduced in toto and also using only every fourth 
data point so that the increment of £ was 10 . Only equation. (17) 
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was. used. The results are presented in table 5«. It can he seen that 
no significant improvement in the variance is obtained by considering 

L , 

the additional data points. 3 Apparently , the distribution of 
crystallites in specimens Wo. 1 and 2 of grade 2BE is not truly symmetri- 
cal about any axis. Specimen No. 2 from the same block of graphite does, 
however, exhibit a symmetry axis. The cause of this behavior is 
probably related to the cause of the variation of the direction of 
the symmetry axis but is not presently known. In any event, it appears 
that .the results for a sample taken at -one location in a block of 
commercial .graphite do not, necessarily, apply at other locations withir 
the block. 

It should be clear from the above discussion that the Ali method 
cannot be expected to give more than a rough approximation of 
crystallite orientation, since this method involves two fairly large 
samples which must, necessarily, be cut from different locations, and 
also since the method depends upon the assumption of an elliptical 
distribution of crystallites, which has been found not to be correct. 

A comparison of Ali measurements of 1(90) /l(0‘) on grades ATJ, 2 d8d., 
and CMB with the values of B 13 for these grades is given in table' 6. 

As expected, the results do not. agree, very well. 

a Note, also, that the values of the various orientation parameters 
are not significantly changed .by considering the additional data, points. 
Comparison of the values in this table with those in table 2(b) 
illustrate the reproducibility of the technique used. 

°The orientation parameters b and B are essentially least- 
squares estimates of the true value of l(90)/l(0). 



TABLE 5 


ORIENTATION PARAMETERS FOR GRADE 2EE USING 
PEAK HEIGHT DATA AND EQUATION ( 17 ) 


Specimen 

number 

At, deg 

B 

M 

o. s , deg , 

0 S ; de g 

d, deg 

a 2 x 10 ' 

1 

2.5 

0.5^ 

2.9 

103 

87 

3 

. 18.6 

1 

10.0 

0.3U 

2.9 

105 

88 

5 ' 

20.7 

2 

2.5 

0.29 

, 3.4 

97 

84 

8 

2.0 

2' 

10.0 

o'. 30 

3.4 

97 

■ 84 

■8 

2.1 

.5 

2.5 

'0.31 

'3.2 

95 

90' 

8 

6.2 

3 ! 

10.0 

0.31 

3.2 

95 

90 

8 

7.4 








TABUS 6 ' 

COMPARISON OP RESULTS OP ALI AND BACON 'METHODS 
(a), Peak Area. Data 


Grade 

1(90.) /l(0), Ali 

B, Bacon a 


o.% 

0.52 

1 

0.96 

' N.S. b 


0.68 

i 

0.66 


(h) Peak Height Data' 



' l(90),/l(0), Ali . 

B, Bacon a 

ATJ 

0.57 

0.50 ; 

2d8d 

0.99 

0.85 

CMB 

0.52 

0.6i§- 


a Average least squares solution to equation (17) 
for all specimens of the grade. 

^No solution. 












CONCLUSIONS AND RECOMMENDATIONS 


The results of this investigation lead to the following conclusion 

1. Most molded, artificial graphites possess an axis about which 
their crystallites are symmetrically oriented. However, there does not 
have to he -a symmetry axis in all such graphites as indicated in this 
investigation by grade 2BE. 

2. The symmetry axis of a molded graphite does not necessarily 
correspond to the forming axis and, in fact, frequently does not. 

5 . The degree of crystallite orientation, the direction of the 
symmetry axis, and, in fact, whether or not a symmetry axis even exists, 
can vary from one location to another within a graphite body. 

It. The equation I =- A cos 3 ^ + B is capable of describing the 
angular distribution of crystallites about the symmetry axis within the 
limits of experimental error. This equation 'is more, precise than the 
polar elliptical equation proposed by Ali, Fitzer, and Ragoss. 

5 . The angle, on a Bacon diffraction film is related to the 
orientation angle, 0, through the equation 

cos 4 = l s cos T6.T5 0 ± mg sin 76.75 0 sin 2 | 

+ \J 1 - Zg 2 - mg 2 sin 76-75 sin 1 

4 , 

where Zg and m s are the direction cosines of the symmetry axis. to 
which (f> is referred. 





6. The .parameters • B, M, Zg, and mg constitute the minimum 
set of parameters necessary to specify completely the orientation of 
crystallites in molded, artificial graphites. 

7 • The height of a diffraction peak is not necessarily proportional 
to the area under the peak'. Consequently, values of the orientation 
parameters "based on peak height data frequently do not agree with the 
values obtained using peak area data. The precision of peak height 
data is generally greater than that of peak area data with the techniques 
employed in this study. 

The following recommendations are made for future work: 

1. The accuracy with which the direction of the symmetry axis can 
"be determined "by the technique " proposed in this thesis should "be 
checked. One possible method would involve a specimen holder with two 
degrees of rotation. First,, a specimen would be aligned as specified 
by Bacon, and the direction of the symmetry axis would be determined as 
described herein. Then the specimen would be. rotated so that the 
symmetry' axis replaced the pressing axis. A second^ X-ray film would be 
exposed and 5 would be determined. If the theoretical method 
presented herein is valid, 5 should be quite small. 

2. A comprehensive study of the factors affecting crystallite 
orientation should be undertaken with the twin goals of ! (-a) understand- 
ing the orientation characteristics of existing graphites and (b) - 

» 4 / 

making possible the manufacture of future graphites ‘w-ith a wide range of 
orientation properties which are consistent throughout the graphite body. 
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3- A reliable measure of l($) should be' determined. 

1)-. The precision of the Bacon method should be increased either 
by using a very fine grain film, such as Kodak Type R or- by using a 

v 

moving counter, rather than film, as the detector. A counter 'detector 
■would eliminate the delay and error- caused by f ilia- processing and 
microdensitometry. An experimental arrangement involving a counter, 
detector which is adaptable to the Bacon method is described in refer- 
ence l4. 



APPENDIX A 


Consider a variable, r, which is an' explicit function of some 
other variable, t, and of u coefficients, Cj 

r = rCc-L, c 2 , . . ., e u , t) (Al) 

If equation (Al) is linear with respect to the coefficients It may be 
rewritten as 

u 

r *1 (A2) 
<3=1 

If v pairs of values of r and t are known they may be 
substituted into- equation (A2) to yield the. family of equations 

u 

r p = 1 = 1 j 2, ..., v (A-3) 

j-1 

The v equations represented by (A3) may be written as the matrix 
equation 


! = AC (A4) 

»here E = .[rjQ, c = [cj, and- A .(A i3 ] = ..(Ajd*)] 

* 

H ^ — b. 1 , equations (A3) or equation. (a 4) can be solved' uniquely 

for tbe values of Cj by Cramer's Rule 9 . However, if v > u the 

determination of the coefficients is overspecified, that is, various^ 

a Provided that |a] f 0 where j A j is the determinant of the matrix ;A. 

> * 
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combinations of u of the equations (A3) can be solved., and each will 
yield a different set of values for the coefficients. No one set' of 
such solutions is preferable to any other set.. However, a statistically 
preferred set of values can be obtained by applying the method of least- 

3. po 

squares . Using this method it can be shown that the best solution 
for the coefficients, c are the elements of the matrix C given by 
the equation 

C = (a^a) “ 1 A T R (A5) 


The coefficients in equation (Al) can- be determined even if the 
equation is nonlinear by linearizing it as its Taylor series expansion, 
truncating all terms higher than first order. 


Ar ~ 



3=1 



Ac 3 


(A6) 


where Ar = r r r' and Ac^ = c j - e^. The superscripts on the 1 

coefficients refer to approximate values about which the -expansion ii 

■ \ 
performed. Superscripts on the dependent variable and its partial 

derivatives indicate evaluation using the approximate . coefficient s . ‘ 

If v .pairs, of .values of r and t are known, equation (a6) 

generates the family of equations 



^he method of least- squares yields the ''best" statistical solution 
that can be inferred from a set of data providing the errors in the data 
are assumed to be normally distributed. 
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AR = AA AC (A8) 

Equation (a 8) is equivalent to equation (A4) and, thus, the 
equivalent least-squares solution is 

AH = (AA T AA)" 1 AA* AR (A9) 

The values of c. are obtained from the elements of AH by means of 
J’ 

the relationship 

e = Ac. + c° (A10) 

J t) J 

Since the higher order terns of the Taylor series expansion of r 

were dropped in linearizing equation (Al) the values of c. obtained 

from equations (A9) and (A10) .are not the exact least-squares values. 

The exact values can be approached as closely as desired, however, by 

repeated solutions of equations (A9) and (A10) using each set of values 

of Cj obtained as an' improved approximation for the next iteration . 

The mathematical procedure just described is essentially that 

outlined in reference; 80. It is the basis of the Fortran subroutine 

FITALL ^ which was used in this study for the solution of equations (17) 

and (19) • .FITALL requires a suitably constructed main program to supply 

SThis statement is correct provided the solution converges properly. 
This will generally be the case if (l) the form of equation (Al) 
correctly describes the data (2) the data are not too badly scattered 
and (5) the first approximations of the parameters -to be determined 
are reasonably good. 



k9 


the necessary input, output,- and control steps, 
a subroutine to supply the values of r? and 


FITALL also requires 
and the matrix 


inversion subroutine MATRXI. 


A listing of the various programs and subprograms used in this 
study is given in appendix B. DRS1 is the main program for equation (17) 
and COSH is the corresponding subroutine.. DRS2 and ELLIPSE are the 
main program and subroutine, respectively, for equation ( 19 ) • AMAX is 
a function subprogram required in normalizing the raw data. 
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PROGRAM DRS1 ( I NPUT , OUTPUT . TAPE5= I NPUT , TAPE6=OUTPuT ) 

DIMENSION DATA ( 144,3 ) ,'P (20 > ,DP(20 ) , RES { 144 ) « C ! (1 44 ) , X I ( 1 44 ) * 

2 PHI ( 1 44 ) . I ( 144 ),RS IG ( 144 ) .ARRAY (403 1 .CNAME (12) ,RESPCT ('! 44 > ,L (4 )>. 

3 I STORE (144) 

REAL I * I MAX , I ZERO , I STORE 

EQUIVALENCE (XI .DATA ( 1 , 1 ) > , ( I .DATA (1.2)). (RSIG.DATA (1.3)) 

EXTERNAL COSN 

LOGICAL TEST 1-.TEST2', TEST3 • 

COMMON /!/ TEST3 
COMMON COSGAM.COSPHI (144) .Cl 
NAMEL I ST/LDATA/-I . XI .NDA-TA 
NAMEL TST/P ARAM/P , NP 

C )=A.P (2 )=B,P (3 )=M,P (4 ).= COS (ALPHA. S) . P (5) =COS ( BETA . S ) 

PRINT 10 

1‘0 FORMAT ( 1 HI ♦ 21HI = A*COS (PH I ) **M + B) 

C READ AND PRINT CASE IDENTIFICATION 

2 READ (5. 100.) CNAME 

100 FORMAT ( 1 2A6 ) 

IF (EOF .5) 998,999 
999 PRINT 101, CNAME 

101 FORMAT ( 1 HI . 12A6//) 

C... ...READ AND PRINT FIRST APPROXIMATION OF COEFFICIENT 
READ ( 5 , PARAM ) 

PRINT 200, tP : (N), N'=l ,NP ) 

200 FORMAT <1X,-*A.(0) = F9.6, 4X-. *8(0) '= *, F9.6,. 4X. *M ( 0 > = *, 

2 F9.6, 4X, *COS (ALPHA, SO (0 ) = *,'F8'.5» 4X, *COS ( BETA.S )■( 0 ) =*’ *♦ 

3 F8 • 55 /V ) ‘ ;• . 

C • « • • » READ AND NORMALIZE DATA — SET WEIGHTING FACTOR, RSIG. TO' UNITY 
READ ('5.LDAT-A) ‘ 

-IMAX=AMAX( I ♦ NDAT A ) 

DO 84 K=1 ,NDATA 
ISTORE (K ) =1 (K ) 

I (0 = 1 (O/lMAx' 

84 RS I G ( K ) = 1 • 

C INITIALIZE PROGRAM PARAMETERS 

VAR=0 

J=1 

' TEST 1=. FALSE. 

TEST-2=. FALSE. 

TEST3=. FALSE. 

L ( 1 ) =NDAT A 
L (2.) =NP 
L (3 ) =200' 

L < 4 ) =0 
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C.-. • ... SOLVE EQUAT I ON ( 1 7 ) FOR 'LEAST— SQUARES VALUES OF COEFFICIENTS 
CALL FIT ALL C.COSN . DAT A , 1 44 , P . DP , 20 , Q , L , ARRAY * 400) 

c TEST SOLUTION AND PRINT APPROPRIATE ERROR MESSAGE IF NECESSARY 

IF (.NOT. TEST3 ) GO TO ISO 

PRINT 190. L . ( D ( N ) . N=1 ,NP) , (DP(N) , N~ 1 . NP ) . COSGAM . COSPH I ( J ) 

190 FORMAT (1-X, ^UNREASONABLE SOLUTION IN PROGRESS.*// IX. *LIST = *, 

2 414// t:x. *P = *, 5F16.6// IX. *DP = *, 5E16.6// IX. 

3 *COS ( GAMMA , S ) = *, E16.6, 4Xi *COS(PHI) = *, E16.6) 

GO TO 2 

ISO TF (L ( 4')-. LE. 30 GO TO 22 
PRINT 21 , L 

21 FORMAT (IX. ^TROUBLE IN FITALL. LIST a *, 414) 

GO TO 2 

22 IF ( L ( 4 ) .LE.,1 ) GO TO 30 
PRINT 23 

23 FORMAT UX, ^CONVERGENCE CRITERION NOT SATISFIED, FOLLOW' I NG RESULT 
2S USE BEST AVAILABLE PARAMETERS.*/' 1 

C. ... .PERFORM FINAL NORMALIZATION 
30 IZERO=P (1 )+P (2 ) 

P-( 1)=P( 1 ) / 1 ZERO 
p (2 )=P ( 2 ) / 1 ZERO 
DO 51 J=1 .NDATA 
I ( J) = I { J)/rZERO 
Cl (J)=CI (J ) /'I ZERO 
RES ( J ) =C I (J)— I ( J) 

RESPCT (J)=RES(J)*100«/I ( J) 

VAR=VAR+RES ( J I **2/ (NqATA-NP ) 

.....TEST COS ( PH I ) FOR UNREASONABLE VALUE AND PRINT VALUE AND ERROR 
....MESSAGE IF FOUND 

‘IF (COSPHI( J). .LE. 1 . 0,)' GO TO 160 
PHI ( JT=COSPHI ( J ) 

TEST1 =. TRUE . 

GO TO 51 

160 PH I { J ) =S7 . 2958*AC0S ( COSPH » ' < ' 1 
51 CONTINUE 

IF (TES'Tl )' PRINT 170' 

170 FORMAT (IX, *ERROR — COS(PHI) GREATER THAN ONE..**) 

....TEST DIRECTION COSINES FOR UNREASONABLE VALUES AND PRINT 
....APPROPRIATE ERROR MESSAGE IF NECESSARY - 
SUMSQ=P ( 4 ) **2+P (5 ) *■*£ 

TF (SUMSO »LE.' 1.0) 'GO TO* 90 
PRINT 110, SUMSQ 

110 FORMAT (IX. 42HERR0R — COS (ALPHA , S4 **2 + COS (BETA , S 1 **2 = 

90 IF (*ABS'(P (4)1 ,GT. 1.0) TEST2= .-TRUE. '* 


i FB.5//) 
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IF ( ABS ( P ( 5 ) ) .GT. I .0 > TEST 2=.. TRUE. 

IF (ABS (COSGAH ). .GT. 1.0) TEST2= . TRUE. 

IF (.NOT. TEST2-) GO TO 120 
PRINT 130 

130 FORMAT (IX. *ERROR — DIRECTION COSINE GREATER THAN ONE.*//) 

GO TO 140 

C. ... .‘CALCULATE ALPHA (S), RETA (3 ) . GAMMA (S). AND DELTA 

120- COSDEL=— COS ( 77 »*1 .74R3E-2 )*P (4 ) + COS ( 13.* 1 .7453E-2 )*COSGAM 
DELTA=57.295S*ACOS(ARS(COSDEL) ) 

IF (COSDEL.LT.O ) DEL TA=1 80 .-DELTA 
ALPHA =57. 2.958* ACOS ( ArS ( P ( 4 ) ) ) 

IF (P(4).LT.O) ALPHA=180. -ALPHA 
BETA =57 . 2958* A COS < ABS ( P ( 5 ) ) ) 

IF (P(5).LT.0) BETA= 180. -BETA 
SGAMMA=57.2958*ACOS (COSGAM) 

C PRINT PROGRAM TEST PARAMETER. LIST 

1’40 PR.INT 20. L 

20 FORMAT (IX. *L IST = *. 4 1 4//) 

C PRINT RESULTS ' - 

PRT NT 40, VAR 

40 FORMAT (IX.- *VARIANCF = *♦ E14.6//) 

PRINT 50, PCI ) ,P.(2) ,P(3) 

50 FORMAT (1X.*A = *. F9.6. 4X , *B = *, F9.6, 4X» *M = *, F9.6//) 
PRINT 60, P(4 ) ,P(5) .COSGAM 

60 FORMAT (IX. *COS(ALPHA,S-) = *» F8.5, 4X. *COS ( BET A , S‘> = *, F8.5, 

2 4X, *COS (GAMMA ,'S ) = *, F8.5//) 

IF { TEST2 ) GO TO 150 

PRINT 70. ALPHA. BETA, SGAMMA, DELTA 

70 FORMAT (IX, *ALPHA(S) = *, F9.4, 4X, *BETA(S) = *, F9.4, 4X. 

2*GAMMA(-S) = *, F9.4, 4X, *DELTA = *, ‘F9.4//1, 

150 PRINT 80, (PHI t-J ) , XI ( J ) ,-I STORE (J),I(J)«CI(J.)« RES ( J ) ,RESPCT ( J ) . 

2 J=1 , NDAT A > 

80 FORMAT -( 37X , *NORMAL I-ZED*/ 5X, *PHI*, 6X, *XI*. BX. *I(EXP>*. 

2 5X, * I (EXP ) * , 4X, * I ( CALC ) * » 4X, *RES* . 4X . *RES(PCT>*/ 

3 (2F10.4, F12.6., 4F10-6)) 

GO TO 2 

998 STOP 
FND 
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SUBROUTINE COSN (RES ♦ DATA 1 NDMAX 1 P « DP <V< NP « J ) 

REAL I 

D I MENS I ON- DATA { NDMAX , 3 )-. P ( 20 ) . DP ( 20 > 

LOGICAL TEST3.TEST4 
COMMON /l/ TEST3 

COMMON COSGAM . COSPH I ( 1 44 ) « C 1(1 44 ) 

TEST4=. FALSE. 

XL I M= 1 .OE-IOO 

C TEST COEFFICIENTS FOR UNREASONABLE VALUE AND SET INDICATOR IF FOUND 

IF <ABS(P(1)> . GT • 2.0) GO TO 20 
IF (ABS(P(2)) .GT. 2.0) Go TO 20 
IF ( ABS ( P ( 3 ) } .GT. 10.0) GO TO 20- 
IF ( ABS ( P { 4 ) ) .GT. 1.0) GO TO 20 
IF ( ABS IP (5 ) ) .GT. 1.0‘) GO TO 20 

C CONVERT XI TO RADIANS 

X I,=OAT A(Jt 1 > 

RX1'=XI *1 . 7453E-2 
I =DATA ( J ,2 ) 

V=DATA-(J,3) 

C CALCULATE COS (BETA. N.) 

ARG1 =S IN (76.75*1 .7453E-2) **2* ( 1 • -SIN (RXI >**2 ) 

IF (ARGI .GT. XL I M ) GO TO’ 60 
COSBETA=0 . 0 
GO TO 70 

60 COSBETA=SQRT ( ARGI ) 

C.....SET SIGN OF COS(BETA-,N) 

if- (-XI .GT. 90. .AND, XI .L7. 270.) C0SBETA=-COSBETA 

CALCULATE COS(GAMMA.S) 

70 ARG2 = ABSt I . -P (4 )**2-p(5)**2 ) 

IF (ARG2 .LT.- XL I M) 60 TO 20 
COSGAM = SORT ( ARGH ) 

IF' . ( C OS GAM .GT. 1.0) GO TO 20 

C CALCULATE COST ALPHA , N 5 

COSALPH=COS ( 76 . 75*1 .7453E-2) 

C. .. ..CALCULATE COSIPHI) 

COSPH I (,J)=COSALPH*P (-4 )+C0SBETA*P<(5 )+COS ( 13.25*1 .'74 53E-2 )*C0SGAM- 
2 *S I N ( RX I ) ■ *’ < 

TEST COST PH I) FOR UNREASONABLE VALUE AND SET INDICATOR IF FOUND 

IF (ABS (COSPHI ( J ) ) .GT. 1.0) GO TO 20 . 

TAKE ABSOLUTE VALUE OF COS (PH I) AND SET INDICATOR FOR- CORRECT' 

SIGN OF DERIVATIVES 
IF ( COSPH I ( J )■ .GE. 0) GO TO 10 
COSPHI ('J ) =— COSPH I (J) 

TEST 4= . TRUE • 
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10 CONTINUE 

...CALCULATE PARTIAL DERIVATIVES OF I (PHI) WITH RESPECT TO COEFF I C I ENTS - AND 
...SET CORRECT SIGN 

DP ( 1) =COSPH I C J ) **P (3 ) 

IF CABS ( DP ( 1 ) ) .LT. xt- I M ) DP(1>=0.0 
DP (2 ) = 1 . 

IF ( ABS ( COSPH I ( J ) ) . GT • XL I M ) GO TO 30 

COSPHI (J)=O.C 
DP (3) =0 . 0 
GO TO 40 

30 DP (3)=P ( 1 ) *COSPH I ( J )**P (3 )*ALOG ( COSPH I ( J } ) 

40 DP (4 )=P( 1 )*P <3)*C0SPHl (J)**(P(3)-1. )* (COSALPH- (P (4 ) /COS GAM)* 

2 COS< 13.25*1 .7453E-2)*SIN(RXI )) 

IF (ABS ( DP ( 4 ) ) .LT. xLI-M) DP(4)=0.0 
IF (TEST4) DP(4)=-CO(4) 

DP (5>=P ( 1 >*P (3 )*C0SPH'I ( J ) ** ( P ( 3 ) — 1 .)* (COSBETA- (P (5 )/C0SGAM->* 

" 2 C0S< 13.25*1 .7453fE-2 )*SIN (RXI > ) 

IF ( ABS ( DP (5)1 .LT. XLIMO DP(5)=0.0 
, IF ( TEST 4 > DP (5 ) =-DP (5 )• 

C CALCULATE KPHIl'IO) 

Cl (J)=P(l }.*C0SPH I <J>**P<3)+P(2> 

C CALCULATE RESIDUAL OF I (PHI) 

RES=C 1 < J ) — I 

c TEST RESIDUAL FOR UNREASONABLE VALUE AND SET -INDICATOR IF FOUND 

IF (ABS (RES) .GT. 1.0) GO TO 20, 

IF ( ABS (RES ) .LT. XLIM) RES=0.0 
GO TO 50 
20 TEST3= . TRUE. 

50 RETURN 
END 
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PROGRAM DRS2 ( INPUT « OUTPUT ,TAPE5= INPUT < TAPE6=OUTPUT ) 

DIMENSION DATA <40.3) ,P<20) .DP (20) . L < 4 ) , RES 1 40 ) . C-I ( 40 ) . X I ( 40 ) . 

2 PHI C40) ..1(40) .RSIG(40 ) . ARRAY <403 ) ,CNAME< 12) . RESPCT (40 ) • 

3 I’ST0RE.(40) 

REAL I . I MAX. I ZERO, (STORE 

EQUIVALENCE (XI .DATA ( 1 , 1 ) ) , ( I , DA TA ( 1' , 2 ) ) . (RSIG.DAT A (1,3)) 

EXTERNAL ELLIPSE 

LOGICAL TEST1 .-TEST2 .TEST3 

COMMON /\/ T.EST3 

COMMON COSGAM.COSPHI (40)-, Cl 

NAMEL I'ST/LDATA/ii ,-XI .NDATA 

NAMELIST/PARAM/P.NP 

C P( 1 )=A,P(2r)=B,P (3 )=COS.( ALPHA ,S') , P ( 4 ) =COS < BET A . S ) 

PRINT 10 

10 FORMAT (TH1, -46HI- = A *8 /SORT ( A**- 2 + (B**2 - A**2 ) *COS ( PH I ) **2 ) ) 

READ AND PRINT CASE IDENTIFICATION 

2 READ (5, 100) CNAME 

100 FORMAT ( I 2A6 ) 

IF (EOF, 5) 998,999 
999 PRINT 101, CNAME 

101 FORMAT (TH1 » 1-2A6//) 

C READ AND PRINT FIRST APPROXIMATION OF COEFFICIENTS 

READ (5.PARAM ) 

P( 1 ) = 1 .00 
NP-4 

PRINT 200, (P(N), N= 1 ♦ NP ) 

200 FORMAT (IX, *A(0) = *, F9.S, 4X. *B(0) = #* F9.6, 4X, 

2 *COS (ALPHA, S-) (O' V = *, F6.5, 4X, *COS (BET A ,S ) ( 0 > = *. F8.5//)_ 
C.....READ AND NORMALIZE DATA — SET WEIGHTING FACTOR, RSIG. TO UNITY 
READ L5.LDA.TA) 

I MAX=.AM AX ( I , NDATA ) 

DO 84 K=l, NDATA 
I STORE (K') = l (K) 

I (K)=I ( K ) / 1 MAX 
84 RS IG (K ) =1 « 

C INITIALIZE PROGRAM PARAMETERS 

VAR- 0 
J=1 . 

TEST 1=. FALSE. ‘ 

TEST2=» FALSE » 

TEST3=. FALSE. 

L < 1 > =NDATA 
L ( 2 ) =NP 
L ( 3) =200' 
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t_.< 4 > =0 

c SOLVE EQUAT I ON (19) FOR LEAST-SQUARES VALUES OF COEFFICIENTS 

CALL F I TALL ( ELL I PSE, DAT A, 40, P, DP, 20, Q,L, ARRAY, 400 ) 

C..... .TEST SOLUTION AND PRINT APPROPRIATE ERROR MESSAGE IF NECESSARY 
IF (.NOT. TEST.3 ) GO TO 180 

PRINT 190, L * (P (N ) , N=1 ,NP ) , (DPIN) , N = 1 ,NP ) , COSGAM , COSPH I ( J ) 

'190 FORMAT (IX, ^UNREASONABLE SOLUTION IN PROGRESS*// IX. *L'I ST = *, 

2 4 14// IX, *P = *, 4E1'6.6// IX, *DP = * , 4E16.6// IX, 

3 *COS ( GAMMA , S ) = *, E16..6, 4X, ’ *COS'(PH I ) = * , E16.6) 

GO TO 2 

180 IF (L ( 4 )-»LE • 3 ) GO TO 22 
- PRINT 21 , L 

21 FORMAT (IX, *TROUBLE IN FITALL. LIST = *, 414) 

* GO TO 2 

22 IF (L(4) . LE . 1 ) GO TO 30 
PRINT 23 

23 FORMAT (IX., *CONVERGENCE CRITERION NOT SATISFIED', FOLLOWING RESULT 
2S USE BEST AVAILABLE PARAMETERS.*//) 

c perform final normalisation 

30 1 ZERO=P ( 1 ) 

P( 1 )=P ( 1 ) / 1 ZERO 
P (2 ) =P ( 2 ) / 1 ZERO ' 

DO 51 J=1 ,N0ATA 
I ( J) = I ( J ) / 1 ZERO 
C‘I (J)=CI ( J ) / I ZERO 
RES (J) = CI'(J> — I ( J) 

RESPCT(J)=RES(J)*100./I ( J) 

VAR=VAR+RES ( J) **2/.( NDATA-NP ) 

....TEST COS (PH I) FOR UNREASONABLE VALUE AND PRINT VALUE AND ERROR ' 
....MESSAGE IF FOUND' 

IF (COSPHI(J) .LE. 1.0) GO TO 160 
PHI ( J)=COSPHI (J) 

TEST 1=. TRUE. 

GO TO 51 

160 PHI <U)=57.29S8*ACOS(COSPHI (J) ) 

51 CONTINUE 

IF (TEST1 ) PRINT 170 

170 FORMAT (IX, *ERROR — COS(PHI) GREATER THAN ONE.*) 

.....TEST DIRECTION COSINES FOR UNREASONABLE VALUES AND PRINT 
.....APPROPRIATE ERROR MESSAGE IF NECESSARY 
SUMSQ=P(3 )**2+P (4 )**? 

IF (SUMSO- .LE, 1.0) GO TO 90 
PRINT 110, 5UMSQ 

no FORMAT (IX, 42HERROR — COS (ALPHA , S ) **2 + COS ( BETA , S ) **2 = 


, F8.5//1- 



90 IF ( ABS ( P ( 3 ) ) .GT. 1.0) TEST2=.TRUE. 

IF {/,'BS CP (4.) ) .GT. 1.0) TEST.'2= .TRUE • 

IF CABSKOSGAM ) ..GT. 1.0) TEST2= .TRUE. 

IF {.NOT. TEST2 ) GO TO 120 
PRINT 1'30 

130 FORMAT (IX. *ERROR — DIRECTION COSINE GREATER' THAN ONE.*//) 

GO TO 140 

....CALCULATE ALPHA (S). BETA(S). CAMMACS). ANO DELTA 
120 COSDEL =— COS ( 77 .*1 .7453E-2 )*P (3 J+COS ( 1 3.* 1 .74S3E-2-) *COSGAM 
0ELTA=57.2958*AC0S ( ABSTCOSOEL) ) 

IF (COSDEL.LT. O') DELTA= ISO. -DELTA 
ALPHA =57 . 2958* AC OS' ( ABS ( P ( 3 ) >') 

IF ( P ( 3 ) .LT. 0) ALPHA=1 80. -ALPHA 
BETA = 57 .2958*AC0S CABS(P(4 ) ) )• 

IF (P ( 4 ) .LT. 0) BETA=180.-BETA'-- 
SGAMMA = 57.2958*A'C0S (COSGAM) 

....PRINT PROGRAM TEST PARAMETER. LIST 
140 PRINT 20, L 

20 FORMAT (lx, *L I ST = *, 4 14//) 

....PRINT RESULTS 
•PRINT 40, VAR 

40 FORMAT (lx, *VAR I ANCE = *, E14.6//) 

PRINT 50, PCI > ,'P(2) 

50 FORMAT ( IX, *A = *, F° .6 , AX', *S = *, F9.6/-/}. 

PRINT 60, PC3).,P<4) .COSGAM 

60 FORMAT (IX, *COS ( ALPHA , s' ) =- *, F8.5, 4X, *COS(BETA,S) = *. F8 
,2 4X ♦ *COS (GAMMA ,S ) = *, F8.5//) 

IF CTEST2') GO TO 150 

PRINT 7.0, ALPHA, BETA, .SGAMMA-, DELTA 

70 FORMAT (IX, *ALPHA(S) = *, F9.4, 4X« *BET A ( S ) = *, F9.4, 4X, 
2*GAMMA(S) = *, F9.4, 4X, *DELTA ='*, F9.4//> 

150 PRINT 80, (PHI (J), XI (J),I STORE (J),I ( J) , C I (J ).,RES ( J ) ,RESPCT { J )‘ 
2 J = 1 , NDATA") 

80 FORMAT (37Xr *NORMAL IZED*/ 5X‘, *PHI*, 6X, *X!*, 8X, *ICEXP)*» 

2 5X , * I (EXP ) * , 4X , * I ( CALC ) * » 4X , *RES* , 4X, *RES(PCT)*/ 

3 (2F10.4, FI 2.6, 4F10»6)) 

GO TO 2 

998 STOP 
END 
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SUBROUT 1 NE ELL I PSE { RES , DATA , NDMAX'.P « DP ,V , NP , J ) 

REAU I 

DIMENSION DATA ( NDMAX , 3 ) , P (.20 ) ,DP ( 20 ) 

LOGICAL TEST3 
COMMON /!/ TEST3 

COMMON COSGAM , COSPH I c40>,CI (40) 

XLIM=1 .OE-IOO 

C.....TEST COEFFICIENTS FOR UNREASONABLE VALUE. AND SET INDICATOR IF FOUND 
IF (ABS(P(l)J .GT. 2.0.) GO TO 20 

IF' ( ABS (P ( 2 ) ) ’.GT. 2.0) GO TO 20 

IF ( ABS ( P ( 3.) ) .GT. 1.0) GO TO 20 

IF ( ABS < P (4 ) ) .GT. 1.0) GO TO 20 

CONVERT XI TO RADIANS 

XI =DAT A { J * 1 > 

RX I =X I *1 .7453E-2 
I =DATA ( J.2 ) 

V=DATA(J,3) 

C. ... .CALCULATE COS (BETA. N) 

ARG1=SIN(76.7S*1 . 7453E-2 )**2* ( 1 .-SIN(RXI )**2) ' 

IF ( ARG1 .LT. XL I M )' GO TO 20 
COSBET.A=SQR.T (ARG1 ) 

C SET SIGN OF COS, (BETA, N) t 

IF (XI .GT. 90. .AND. XI .LT. 270.) COSBETA=-'COSBETA 
C...... CALCULATE COS (GAMMA, S)’ 

ARG2=ABS ( 1 ,-P(3)**2-P(4)**2) 

IF (-ARG2 .LT. XLIM) GO TO 20 
COSGAM=SQRT ( ARG2 ) 

IF (COSGAM .GT • 1.0) GO’ TO 20 

CALCULATE COS (ALPHA, N) 

•C0SALPH=C0S'(76.75’*1 .74S3E-2) ' 

CALCULATE COS (PHI) 

COSPH I { J)=C0SALPH*P(3 >+COSBETA*P (4 >+COS’( 1 3.25#-l .7453E-2 )*C0SGAM’ 

2 *S I N( RXT ) 

C . . . . . TEST COS (PH I) FOR UNREASONABLE VALUE AND SET INDICATOR IF FOUND 
IF (ABS (COSPHI ( J) ) .GT. 1.0) GO TO 20 

CALCULATE INTERMEDIATE FUNCT TONS TO BE USED BELOW 

S I NSQ= ( 1 .-COSPH 1 ( J ) **2 ) 

ARG3=P-(1 )**2+(P(2)**2-P<l )**2 ) *COSPH I (J)**2 
IF’ ( ARG3 .LT. 0.) GO TO 20 
VAL1=SQRT (ARG3 > 

IF (VAL1 . LT . XLIM) GO TO 20’ • 

VAL2=-P ( 1 )*P(2)*(P<2)**2-P( 1 )**2)*C0SPHI ( J) 

VAL3=C0S-( 13.25*1 .7453E-2 ) *S I N ( RX I )/COSGAM 
C. ....CALCULATE PARTIAL DERIVATIVES OF I(PHI) WITH RESPECT TO COEFFICIENTS 
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DPU )=PC2)/VAL1-P-(1 )**2*P(2>*SINSQ/VALI**3 
IF CABS CDP ( I ) ) »LT . XL'IM) DP(I)=0_.0 
DP C2‘) =P Cl )/VALl-P(l )*P (2 }**2*C0SPHI ( J")**2/VAUI **3 
IF CABSCDPC2)) .LT. XLIM) DP(2)=0.0 
DP(3-)=VAL2/VAL1**3*CC0SALPH-VAL3*PM ) } 

IF ( ABS ( DP ( 3 > ) »LT. XL I M ) DP(3)=0.0 
DP ( A ) =VAL2/VAL 1 **-3* C COSBETA~VAL3#P (2 > } 

IF (ABSCDPC4)) »LT. XL I MO DP(4) = 0.0 
C. ....CALCULATE I (PHI) CO) 

Cl <J)=PCi )*PC2)XVALI 
C. ... .CALCULATE RESIDUAL OF I (PHI) 

RES=C I C.J ) — I 

C..L..TEST RESIDUAL FOR UNREASONABLE VALUE AND’ SET INDICATOR IF FOUND 
IF (ABS(RES) .GT.'l.-rv) GO TO 20 
IF (ABSCRES) .LT. XLIM) RES=0.0 
GO TO 50 
20 TEST3=.TRUE. 

50 RETURN 
END 
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SUBROUTINE F I TALL (RES i 0 4 X , NDMAX,. p,DP .MAX .Q .L I ST 4 HOLD .MMAX ) ' 
DIMENSION PC 20) 4 DP (20 ) »LIST.(4)-.A (20»20) .S C400) <H0LD.C403 ) .B’C-20 ) « 
.2 ISO) . ' ' * , 

EQUIVALENCE ( IS ( 1 ) ,Kk) . ( IS('2 ).»NF IX) . C IS (31 , KS IG ) . < S C 1 ) ♦ A C 1 . 1 ) ) 
DATA CONV/O • 00 1 / 

LOGICAL TEST3 
COMMON /l/ TEST3 
C..... INITIALIZE 
JJ=L I ST <!• )•» 

KK=LIST(2). 

LOOPS=L I ST ( 3 ) 4 
L I ST f 3 ) =0 . 

NFREE = JJ-KK 4 

IF (LOOPS.EQ.O )G0 TO IO 4 

MOST=2*KK4 

IF (KK.LT. 1 .OR.NFREE.'lT.OIGO TO 14 
!FCKK.GT4MAX)GO TO 154 
C BEGIN ITERATION LOOP 

1 nfix=ki<4 

DP' 2 N=1 , MAX 4 

B(N)=0.4 

DO 2 M= 1 4 MAX 4 

2 A(H,N)=0.-, 

DO 6 J = 1 4 JJ 4 

CALL RES I D ( R 4 X 4 NDM AX 4 P 4 DP 4 VAR 4 KIC. 

'IF (TEST3 ) GO TO 29 
IF ( VAR 4 NE 4 0’4 ) GO TO 4, 

NFIX=NFIX+1 4 
IFCNFIX.GT-.MAXJGO TO 154 

ifcnfix4GT4MOst>go to 144 

c. 4 4 . 4 CALCULATE FIX POINT FLEMENTS 
.BCNFIX)=R 
DO 3 M=l 4KK 
A C M 4 NF I X ) =DP CM ) 

3 ACNFIX4.M).= DP;CM) 

GO TO 6 

C4 4 . 4 .CALCULATE REGULAR POINT 

4 DO 5 N=1 .KK 

B CN) =B (N)+R*DP CN)/VAR4 
DO 5 M=N 4 KK 

5. ACM,N)=A(M4N)+DPCM)*DPCN l/VAR. 

6 CONTINUE 

C FINISH OFF MATRIX 

DO 7 N=1 ,KK 
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DO- 7 M = N,KK 

7 A(N.M)=A’(M*N> . 

-INVERT MATRIX 

KSIG=0, 

CALL I MAT INV t A * MAX* NF TX >DETi KS I G ) 
IFtKSIG.EQ.l *OR.KSIG.EQ.2)GO TO ‘l6, 
IFCLOOPS.EQ.-l )GO TO 1<0. 

CALCULATE NEW VALUES OF PARAMETERS 

B I G=0 • * 

DO 9 M=1 iKK 

SUM= 0 • i 

DO a N— 1 * NF I’X 

8 SUM=SUM— A"{ M'» N ) *8 ( N ) * 

PCM)=P(M.)+SUM» 

9 BIG=9IG+AMAX1 (ABS(SUM)-CONV*ABS(P<M) )*0. >* 
LISTC3 )=LIST13> + 1 t, 

IF (BIG.EO.'O. )GO To' 10*. 

I F ( L I ST C 3 ) • GE • LOOPS ) GO TO 17. 

GO TO 1 

C. . . ...CONVERGENCE ACHIEVED — COMPUTE VAR I 

10 IF (LOOPS. EG. 1 )GO TO 18 * 

Q = 0. « 

SUM=0. * 

DO 11 J=1*JJ. 

CALL REST D(R*Xi NDMAX * P * DP *.VAR*KK « J ) 

IF ( VAR.EO.O. >GO TO ll'« 

SUM=SUM+R*R/VAR * 

11 CONTINUE 

Q=SQRT (SUM/FLOAT CNFREET ) * 

IFCKSTG.NE.OIGO TO 19* 

GO TO 50* 

14 L I ST'C 4 ) =L I ST C4 ) +3 , 

15 LIST (4 )‘=LIST (4 > + l , 

16 L-I ST ( 4 ) =L I ST ( 4 ) + 1 * 

17 L I ST.t 4 ) =L I ST •( 4 ) + 1 * 

18 LIST (4 >=LIST t.4 ) + l * 

19 L 1 ST (4 )=L 1ST (4 1 + 1 * 

C ENTRY FOR SAVING MATRl-X 

ENTRY FITSAV 
50 DO 20 N=1 *MMA'X 

20 HOLD {N+3 ) =S ( N-) 

DO 27 N= 1 *3 

27 HOLD (N)=ISCN) 



...ENTRY FOR RESTORING SAVED MATRIX 
ENTRY FITBAK 
DO 21 N= 1 » MMAX 
21 SCN)=H0LD(N+3) 

DO 2B N=1 ,3 
26 IS(N)=HOLD(N> 

29 RETURN 
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SUBROUTINE MATRXI < X-,MMAX ( MX . DET 4 KS I G ) • 

C SUBROUTINE FOR OBTAINING DETERMINANT OF INPUT MATRIX 

D I MENS I ON, X (-20 4 20.) , Y ( 20 , 20-) 4 Z < 2 »,20 . 20 ) 4 K ( 20 ) 

LOGICAL MATIN 
EQUIVALENCE (Y,Z) 

DOUBLE PRECISION A, R, PROD, .Y 4 

DATA DELT 4 EPS 4 LOOPS/0 • 000 1 4 1. .E-~, », 

GO TO 71 1 

MATRIX INVERSION BY GAUSS-JORDAN ELIMINATION 

ENTRY IMATI'NV 
MAT I N=. TRUE. 

GO TO 1 

711 MAT I N=. FALSE-. 

1 DO 405 1=1.2 
DO 405 J= I 4 20 

DO 405 I CNT = 1 420 
Z ( I 4 J. 1 CNT ) =0 • 

405 CONTINUE 

INITIALIZE ROUTINE AND TEST MX < =ORDER OF MATRIX) 

M=MX. 

.IF (M.GT . I .AND.M.LE.MMAX )G0 TO 5. 

IFtM.EQ.l ) GO TO 2. 

KSIG=KSIG+1 4 
RETURN 

2 PROD=X ( 1 4 1) 4 

I F (PROD.NE. 0 . ) GO TO 44 

3 KS I G=KS IG+24 
RETURN 

4 X< 1 4 1 >=1 ./PROD 1 . 

GO TO 23 

5 PR0D=1 . . 

MM=M— 1 4 - 

DO 6 1 = 1 4 M> 

K(I) = I 4 

DO 6 J=1 <M< 

6 Y ( I 4 J)=X( I-. J) 4 

C.... -.BEGIN BY FINDING LARGEST PIVOTAL ELEMENT 
DO IT 1=1 4 M 4 
A = 0 . . 

DO 7 J=I 4 M 4 

I F CDABS ( Y ( J , 1 ) ) .LE. A )GO TO 7 
A=DABS (Y < J. 1 ) ) 

L=J. 

7 CONTINUE 
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IF (A.EQ.'O. > GO TO 3* 

C. . .-. .REARRANGE ROWS AND ORDER ARRAY 
N=K(L), 

KCL)=K( I ) , 

KM ) =N * 

DO 8 J=1 ,M, 

A=Y( I , J), 

Y( I . J)=Y(L»J ) . 

8 Y(L,J)=A, 

C..... REDUCE PIVOTAL ROW 
A = Y (1,1). 

I F ( .NOT . MATT N ) PROD=PROD*A 
DO 9 J= 1 .MM. 

9 Y( I . J)=Y( I , J+l )/A. 

Y ( 1 ♦ M ) = 1 ./A. 

REDUCE REMAINING ROWS 

DO 11 L=I,M, 

IF (L.EQ. I )GO TO 11, 

A=YCL, 1 ) ♦ 

DO 10 N= 1 ,'MM , 

Y (L.N )=Y (L.N+l )-A*Y< I , N ) , 

IF ( DABS ( Y (L » N ) ) • L T • { DABS ( Y (L i N+l ) ) *EPS > )Y(L.N)=0. 
. 10 CONTINUE, 

Y (L,M >=-A*Y ( I ,M ) ♦ 

11 CONTINUE 

C... .. .UNSCRAMBLE INVERTED MATRIX 
DO 15 1=1, M, 

I F (K ( I ) .EQ. I .) GO TO 15. 

PROD=— PROD., 

DO 12 J=I,M, 

IF(K(J).EQ.I ) GO TO 1.3, 

12 CONTINUE. 

GO TO 3 

13 DO I A L = 1 ,M, 

A = Y(L, I 1 , 

Y ( L , I )=Y(L,J) , 

14 Y(L,J)=A, 

K(J)=K(.I > , 

15 CONTINUE 

C OBTAIN ERROR MATRIX 

DO 20 N=-l, LOOPS, 

TEST = 0. , 

DO 17 1=1 ,M, 

DO 17 J=1 ,M, 
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R= 0 

DO 16 L = 1 .M, 

16 R=R— Z <1 iLi J! (I <L) i 
IF (1 .EQ. J’)R = R+1 • * 
ABSR=SNGL (DABS (R ) ) 

TEST = AMAX1 ( TEST « ABSR ) 

17 Z ( 2 « NJ)=Ri 
DO 1-9 1=1'* M, 

DO 19 J=1 iMi 
A=0. * 

DO 18 L=1 iMi 

18 A = A+“Z (1 ( I (L) *1 (2 .L . J ) » 
ZU(l(J) s Z(lil(Jl+Ai 

19 CONTINUE. 

I F.( TEST »LE* DELT 1GO TO 21. 

20 CONTINUE 

KS I G=KS I G+3 » 

transfer FINAL inverse 

21 DO 22 I =1 .M, 

DO 22 J=l.Mi 

22 X < t . J ) =Z ( 1 . I . J ) . 

23 DET=PROD . 

RETURN 

END 
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FUNCTION AMAX ( A »N ) 
DIMENSION A (144) 
AMAX=A( 1 ) 

DO 70 M=2,N 
IF (A (M > -AMAX >70 .70,60 
60 AMAX=A(M) 

70 CONTINUE 
RETURN 
END 
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